{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f3e24f2e",
   "metadata": {},
   "source": [
    "# Randomized Benchmarking\n",
    "\n",
    "\n",
    "*Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved.*\n",
    "\n",
    "\n",
    "## Outline\n",
    "\n",
    "**Note: The number of points in the Quantum Hub account and time spent to run this tutorial program will vary based on the parameters users input. Users usually need about half-hour and 100 points to obtain relatively reliable results. If you want to get more points, please contact us on [Quantum Hub](https://quantum-hub.baidu.com). First, you should log into [Quantum Hub](https://quantum-hub.baidu.com), then enter the \"Feedback\" page, choose \"Get Credit Point\", and input the necessary information. Submit your feedback and wait for a reply.**\n",
    "\n",
    "There are generally two ways in experiments for characterizing the performance of a quantum computer in a superconducting platform: Quantum Process Tomography(QPT) and Randomized Benchmarking(RB) \\[1\\]. QPT can completely characterize a gate, decomposing a process into Pauli or Kraus operators, but improving gates by QPT is complicated and resource-consuming. Also, State Preparation And Measurement (SPAM) errors can be confused with process errors. However, RB is a concept of using randomization methods for benchmarking quantum gates. It is a scalable and SPAM robust method for benchmarking the full set of gates by a single parameter using randomization techniques. So it is useful to use RB for a relatively simple and at least SPAM-robust benchmark, especially when the number of qubit increases.\n",
    "\n",
    "In this tutorial, we will implement RB on one of the qubits in our noisy simulator defined in advance to characterize the average error rate on a Hadamard gate. The outline is as follow:\n",
    "\n",
    "- Introduction\n",
    "- Preparation\n",
    "- Define the hardware to benchmark\n",
    "- Running RB\n",
    "- Summary"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "884577b5",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "\n",
    "**Basic RB**\n",
    "\n",
    "We usually use pure Clifford gates to get a sequence fidelity decay curve by randomly choose $m$ Clifford gates followed by an inverse Clifford gate to make sure that the effect of all $m+1$ Clifford gates is an identity operation $\\mathcal{I}$ and then apply them to a qubit:\n",
    "\n",
    "![basicRB](figures/basicRB.png)\n",
    "\n",
    "The above circuit shows how our basic randomized benchmarking works. $C_{i}$ represents the $i_{th} (i = 1, 2, 3, \\dots, m)$ Clifford gates we randomly choose. Ideally, assume our qubit's initial state is $|\\psi\\rangle$, and there is no noise in this RB procedure, we will get the final state being $|\\psi\\rangle$ for $100\\%$ probability, and we use the probability of the final state still being $|\\psi\\rangle$ as a measure of our RB sequence fidelity \\[2\\]. However, in reality, this sequence fidelity $\\mathcal{F}$ decays when the number of gates increases because of the accumulation of noise. If the noise distribution is time- and gate-independent, such decay can be described by a zeroth fitting function:\n",
    "\n",
    "$$\n",
    "\\mathcal{F}^{(0)}=Ap_{\\rm basic}^m+B,\n",
    "$$\n",
    "\n",
    "where $m$ is the number of Cliffords we apply. For more detailed basic knowledge and related theories about randomized benchmarking, users can refer to \\[3\\].\n",
    "\n",
    "\n",
    "The advantage of RB is that the SPAM error is included in the parameters $A$ and $B$ without affecting the rate of decay parameter $p$. More specifically, if we decompose the initial density operator $\\rho$ and measurement operator $\\hat{E}$ into Pauli basis $\\hat{P}_i$: \n",
    "\n",
    "$$\n",
    "\\rho=\\sum_jx_j\\hat{P}_i/d,\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\hat{E}=\\sum_j\\tilde{e}_j\\hat{P}_j,\n",
    "$$\n",
    "\n",
    "then $A = \\sum_{j\\neq 0}\\tilde{e}_jx_j$, and $B = \\tilde{e}_0$ where $d\\equiv{2^n}$ and $n$ is the number of qubits. Once we successfully get $p_{basic}$ by fitting the curve, we can obtain the average EPC(Error-rate Per Clifford) immediately:\n",
    "\n",
    "$$\n",
    "{\\rm EPC}=\\frac{(1-p_{\\rm basic})(d-1)}{d}.\n",
    "$$\n",
    "\n",
    "\n",
    "**Interleaved RB**\n",
    "\n",
    "Interleaved randomized benchmarking is used to get the average error rate of a particular quantum gate. After we obtained a reference sequence fidelity decay curve by first implementing the basic RB, we can then apply the so-called interleaved randomized benchmarking \\[2\\] to benchmark the performance of the Hadamard gate. We interleaved our aimed gate behind every time we implement a random Clifford gate and then make the gate sequence equals to the identity gate by an inverse Clifford gate. The interleaved randomized benchmarking sequence shown below takes the Hadamard gate(H gate) being the target gate as an example: :\n",
    "\n",
    "![interleavedRB](figures/interleavedRB.png)\n",
    "\n",
    "And we can as well get a sequence fidelity decay curve that is similarily described by:\n",
    "\n",
    "$$\n",
    "\\mathcal{F}^{(0)\\prime}=A^{\\prime}p_{\\rm gate}^m+B^{\\prime}.\n",
    "$$\n",
    "\n",
    "Consequently, we calculate the average error rate on the aimed gate EPG(Error-rate Per Gate) by subtracting the reference curve away:\n",
    "\n",
    "$$\n",
    "r_{\\rm gate}=\\frac{(1-p_{\\rm gate}/p_{\\rm ref})(d-1)}{d}.\n",
    "$$\n",
    "\n",
    "As we can see, this gate error $r$ characterizes the gate's performance we want to benchmark since it is derived from a set of experiments' data.\n",
    "\n",
    "Next, we will use Quanlse to implement the RB experiment to benchmark a Hadamard gate. As mentioned above, since we are going to benchmark a specific gate, we need to use both basic RB and interleaved RB. It is worth mentioning that in the up-to-date version of Quanlse, we only support the single-qubit Randomized Benchmarking with Clifford basis. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a287f60",
   "metadata": {},
   "source": [
    "## Preparation\n",
    "\n",
    "First, we need to import the following necessary packages and get a token in order to use the cloud service:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d82829f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the necessary packages\n",
    "from Quanlse.Utils.RandomizedBenchmarking import RB\n",
    "from Quanlse.Utils.Functions import basis, tensor\n",
    "from Quanlse.QOperation import FixedGate\n",
    "from Quanlse.Superconduct.Simulator import PulseModel\n",
    "from Quanlse.Superconduct.SchedulerSupport import SchedulerSuperconduct\n",
    "from Quanlse.Superconduct.SchedulerSupport.GeneratorRBPulse import SingleQubitCliffordPulseGenerator\n",
    "from Quanlse import Define\n",
    "\n",
    "from math import pi\n",
    "from scipy.optimize import curve_fit\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# Import Define class and set the token\n",
    "# Please visit http://quantum-hub.baidu.com\n",
    "from Quanlse import Define\n",
    "Define.hubToken = ''"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8542195",
   "metadata": {},
   "source": [
    "## Define the Hardware to Benchmark"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce848df3",
   "metadata": {},
   "source": [
    "Then we need to construct a virtual quantum hardware - a noisy simulator and decide which of the qubit and gates we want to benchmark.\n",
    "\n",
    "Quanlse supports user-defined multi-qubit noisy simulator, for more detail, users can refer to [multi-qubit noisy simulator](https://quanlse.baidu.com/#/doc/tutorial-multi-qubit-noisy-simulator). Here, we build up a two-qubit system using Quanlse. Each qubit takes three energy levels into consideration. Then, we want to benchmark the performance of the Hadamard gate implementation on the first qubit of this virtual hardware: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa667d7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the basic parameters of the simulator\n",
    "sysLevel = 3  # The number of energy levels of each qubit\n",
    "qubitNum = 2  # The number of qubits simulator has\n",
    "\n",
    "# Qubit frequency & anharmonicity\n",
    "wq0 = 5.282 * (2 * pi)  # The frequency for qubit 0, in 2 pi GHz\n",
    "wq1 = 5.248 * (2 * pi)  # The frequency for qubit 1, in 2 pi GHz\n",
    "anharm0 = - 0.33081 * (2 * pi)  # The anharmonicity for qubit 0, in 2 pi GHz\n",
    "anharm1 = - 0.33043 * (2 * pi)  # The anharmonicity for qubit 1, in 2 pi GHz\n",
    "qubitFreq = {0: wq0, 1: wq1}\n",
    "qubitAnharm = {0: anharm0, 1: anharm1}\n",
    "\n",
    "# Coupling map between qubits\n",
    "g01 = 0.002 * (2 * pi)  # The coupling strength of the interaction between qubit 0 and qubit 1, in 2 pi GHz\n",
    "couplingMap = {(0, 1): g01}\n",
    "\n",
    "# Taking T1 & T2 dissipation into consideration, in the unit of nanosecond\n",
    "t1List = {0: 50310, 1: 62200}\n",
    "t2List = {0: 13630, 1: 24280}\n",
    "\n",
    "# Sampling time\n",
    "dt = 1.  \n",
    "\n",
    "# Build a virtual QPU\n",
    "model = PulseModel(subSysNum=qubitNum,\n",
    "                   sysLevel=sysLevel,\n",
    "                   couplingMap=couplingMap,\n",
    "                   qubitFreq=qubitFreq,\n",
    "                   dt=dt,\n",
    "                   qubitAnharm=qubitAnharm,\n",
    "                   T1=t1List, T2=t2List,\n",
    "                   ampSigma=0.0001)\n",
    "ham = model.createQHamiltonian()\n",
    "\n",
    "# The initial state of this simulator\n",
    "initialState = tensor(basis(3, 0), basis(3, 0))\n",
    "\n",
    "# Decide the qubit we want to benchmark\n",
    "targetQubitNum = 0\n",
    "hamTarget = ham.subSystem(targetQubitNum)\n",
    "\n",
    "# Decide one specific gate we want to benchmark\n",
    "targetGate = FixedGate.H"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a349539",
   "metadata": {},
   "source": [
    "And now we are ready to execute the RB experiment on our chosen qubit.\n",
    "\n",
    "Since we will obtain plenty of pulses when appling RB experiment, we have to use `SchedulerSuperconduct()` to efficiently schedule those pulses:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7458f22",
   "metadata": {},
   "outputs": [],
   "source": [
    "sche = SchedulerSuperconduct(dt=dt, ham=hamTarget, generator=SingleQubitCliffordPulseGenerator(hamTarget))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d321995",
   "metadata": {},
   "source": [
    "## Running RB"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00d80a50",
   "metadata": {},
   "source": [
    "We can obtain a reference curve by doing basic RB calling ``RB`` function which takes some input into consideration: the system ``model`` we have and its initial state ``initialState``, the index ``targetQubitNum`` of the qubit in this system we want to benchmark, the list of Clifford gate number `size` and the number of sequences ``width`` every element in ``size`` has, the scheduler `sche` we initiate before, the sampling time `dt`. We also have to decide whether to use the interleaved RB method. If `interleaved=True`, which means we now implement interleaved RB, then the parameter `targetGate` representing the gate we want to benchmark must be included. In this tutorial, we take qubit's $T_1$ and $T_2$ into consideration, which means that we simulate the RB experiments within an open system, so we set `isOpen=True`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47bccc9d",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Create a list to store the outcome\n",
    "sizeSequenceFidelityBasic = []\n",
    "sizeSequenceFidelityInterleaved = []\n",
    "\n",
    "# Core parameters of an RB experiment\n",
    "size = [1, 10, 20, 50, 75, 100, 125, 150, 175, 200]\n",
    "width = 5\n",
    "\n",
    "# Start RB experiment. First get a basicRB curve used for reference. Then implement the interleavedRB to benchmark our Hadamard gate\n",
    "for i in size:\n",
    "    widthSequenceFidelityBasic = RB(model=model, targetQubitNum=targetQubitNum, initialState=initialState, size=i, width=width, sche=sche,\n",
    "                                    dt=dt, interleaved=False, isOpen=False)\n",
    "    sizeSequenceFidelityBasic.append(widthSequenceFidelityBasic)\n",
    "print(sizeSequenceFidelityBasic)\n",
    "    \n",
    "for j in size:\n",
    "    widthSequenceFidelityInterleaved = RB(model=model, targetQubitNum=targetQubitNum, initialState=initialState, size=j, width=width,\n",
    "                                          targetGate=targetGate, sche=sche, dt=dt, interleaved=True, isOpen=False)\n",
    "    sizeSequenceFidelityInterleaved.append(widthSequenceFidelityInterleaved)\n",
    "print(sizeSequenceFidelityInterleaved)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db1668e7",
   "metadata": {},
   "source": [
    "After we successfully get plenty of sequence fidelity data that these two RB method produced, we then fit the curve to obtain both the average EPC and the average EPG:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02a4a7be",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the fitting function\n",
    "def fit(x, a, p, b):\n",
    "    \"\"\"\n",
    "    Define the fitting curve\n",
    "    \"\"\"\n",
    "    return a * (p ** x) + b\n",
    "\n",
    "# Define the function of calculating the EPG(Error-rate Per Gate) with p_{gate} and p_{ref}\n",
    "def targetGateErrorRate(pGate, pRef, dimension):\n",
    "    \"\"\"\n",
    "    Calculate the specific gate error rate\n",
    "    \"\"\"\n",
    "    return ((1 - (pGate / pRef)) * (dimension - 1)) / dimension\n",
    "\n",
    "\n",
    "# Get the EPC(Error-rate Per Clifford) and p_{ref}\n",
    "fitparaBasic, fitcovBasic = curve_fit(fit, size, sizeSequenceFidelityBasic, p0=[0.5, 1, 0.5], maxfev=500000,\n",
    "                                      bounds=[0, 1])\n",
    "pfitBasic = fitparaBasic[1]\n",
    "rClifford = (1 - pfitBasic) / 2\n",
    "print('EPC =', rClifford)\n",
    "\n",
    "# Get the parameter p_{gate}\n",
    "fitparaInterleaved, fitcovInterleaved = curve_fit(fit, size, sizeSequenceFidelityInterleaved,\n",
    "                                                  p0=[fitparaBasic[0], 1, fitparaBasic[2]], maxfev=500000,\n",
    "                                                  bounds=[0, 1])\n",
    "pfitInterleaved = fitparaInterleaved[1]\n",
    "yfitBasic = fitparaBasic[0] * (pfitBasic ** size) + fitparaBasic[2]\n",
    "yfitInterleaved = fitparaInterleaved[0] * (pfitInterleaved ** size) + fitparaInterleaved[2]\n",
    "EPG = targetGateErrorRate(pfitInterleaved, pfitBasic, dimension=2)\n",
    "print('EPG =', EPG)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "459e9f78",
   "metadata": {},
   "source": [
    "And plot the curve to obtain the result of the decay curve of sequence fidelity for our RB experiment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8279d3fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the decay curve of our RB experiment\n",
    "plt.figure(figsize=(18, 6), dpi=80)\n",
    "plt.figure(1)\n",
    "ax1 = plt.subplot(121)\n",
    "ax1.plot(size, sizeSequenceFidelityBasic, '.b', label='experiment simulation data')\n",
    "ax1.plot(size, yfitBasic, 'r', label='fitting curve')\n",
    "plt.xlabel('$m$')\n",
    "plt.ylabel('Sequence Fidelity')\n",
    "plt.title('basic RB using Quanlse')\n",
    "plt.legend()\n",
    "ax2 = plt.subplot(122)\n",
    "ax2.plot(size, sizeSequenceFidelityInterleaved, '.b', label='experiment simulation data')\n",
    "ax2.plot(size, yfitInterleaved, 'r', label='fitting curve')\n",
    "plt.xlabel('$m$')\n",
    "plt.ylabel('Sequence Fidelity')\n",
    "plt.title('interleaved RB using Quanlse')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e61bd138",
   "metadata": {},
   "source": [
    "Here, $m$ represents the number of Clifford gates we applied. The curve shown in the figure reflects the phenomenon that the accumulated noise as the number of gates (the number of pulses) increases causes the sequence fidelity to decay exponentially. It can be seen that with this scheme, we can automatically generate high-precision pulses that are adapted to the gate operation of the target quantum hardware, perform pulse scheduling when the number of pulses increases significantly, and further conduct randomized benchmarking experiments on the quantum hardware."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7d9047f",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "This tutorial describes how to implement RB experiments to benchmark the performance of a specific gate on one of the qubits in our noisy simulator. Users can click on this link [tutorial-randomized-benchmarking.ipynb](https://github.com/baidu/Quanlse/blob/main/Tutorial/EN/tutorial-randomized-benchmarking.ipynb) to jump to the corresponding GitHub page for this Jupyter Notebook documentation to get the relevant code, try the different parameter values for better curve fitting, and further exploring the variant method of the Randomized Benchmarking."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "214a647a",
   "metadata": {},
   "source": [
    "## Reference\n",
    "\\[1\\] [Kelly, Julian, et al. \"Optimal quantum control using randomized benchmarking.\" *Physical review letters* 112.24 (2014): 240504.](https://doi.org/10.1103/PhysRevLett.112.240504)\n",
    "\n",
    "\\[2\\] [Magesan, Easwar, et al. \"Efficient measurement of quantum gate error by interleaved randomized benchmarking.\" *Physical review letters* 109.8 (2012): 080505.](https://doi.org/10.1103/PhysRevLett.109.080505)\n",
    "\n",
    "\\[3\\][Magesan, Easwar, Jay M. Gambetta, and Joseph Emerson. \"Scalable and robust randomized benchmarking of quantum processes.\" *Physical review letters* 106.18 (2011): 180504.](https://doi.org/10.1103/PhysRevLett.106.180504)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
